      subroutine uprefst(istepp)

      use alloc
      use refstate

      implicit real*8(a-h,o-z)
      character*80 line,fndebug
      integer year1,month1,day1,hour1,min1,sec1
      integer year2,month2,day2,hour2,min2,sec2
      real(kind(0.d0)),allocatable,dimension(:,:,:)::
     :  xflop,yflop,ubxpro1,vbypro1,rhox,rhoy
     :  ,ubxpro2,vbypro2,dubxprodt,dvbyprodt
     :  ,dubxprodz1,dubxprodz2,dvbyprodz1,dvbyprodz2
     :  ,dubxprodzdt,dvbyprodzdt

      if(verbose.ge.2)
     :   write(*,*) 'uprefst:',istepp,intrefst

      if(istepp.eq.intrefst) then
        if(verbose.ge.1) write(*,*) 'uprefst: cycle',istepp

        read(35,'(a)',end=8888) line
        read(35,*) tprof2
        do ip=1,nprof
          read(35,'(a)',end=8888) line
          if(verbose.ge.2) write(*,*) 'Read prof:',ip,line
          read(35,*) xprof(ip),yprof(ip)
          read(35,'(a)',end=8888) line
          do ii=0,ndat-1
            read(35,*) zthdat(ii,ip),thdat(ii,ip)
     :        ,usdat(ii,ip),vsdat(ii,ip),qvsdat(ii,ip),pressdat(ii,ip)
          enddo
          zthdat(ndat,ip)=max(20000.,zthdat(ndat-1,ip))
          usdat(ndat,ip)=usdat(ndat-1,ip)
          vsdat(ndat,ip)=vsdat(ndat-1,ip)
          qvsdat(ndat,ip)=qvsdat(ndat-1,ip)
          thdat(ndat,ip)=thdat(ndat-1,ip)
          pressdat(ndat,ip)=100.  ! not relevant?
        enddo

        zusdat=zthdat
        zvsdat=zthdat
        zqvsdat=zthdat

        fndebug='NewWeipro.grd'
        call weiprof(fndebug)
c
c mrefst=timesteps until next update
c intrefst=nstep of update
c
        if(realtime) then
c         write(*,*) 'uprefst:',tprof1,tprof2
          write(line,'(f15.0)') tprof1
          write(*,*) line
          read(line,'(i4,5i2)') year1,month1,day1,hour1,min1,sec1
          write(line,'(f15.0)') tprof2
          read(line,'(i4,5i2)') year2,month2,day2,hour2,min2,sec2
          h1=julian(year1,month1,day1)
          h2=julian(year2,month2,day2)
          stepseconds=((h2*24+hour2-h1*24-hour1)*3600.
     :      +(min2*60+sec2-min1*60-sec1))
          mrefst=nint(stepseconds/dt)
        else
          mrefst=nint((tprof2-tprof1)/dt)
        endif

        intrefst=istepp+mrefst
        if(verbose.ge.1) write(*,*) 'uprefst: intrefst',intrefst
        write(nchn1,'(a20,2f15.0,2i8)')
     :    'Next bc profile:',tprof1,tprof2,istepp,intrefst

c save current ref profiles

        thpro1=thpro0
        uspro1=uspro0
        vspro1=vspro0
        qvspro1=qvspro0
        press1=press0
        rhopro1=rhopro

        dthdz1=dthdz
        dusdz1=dusdz
        dvsdz1=dvsdz
        dqvsdz1=dqvsdz

        call refprofs

        do ip=1,nprof
          dpressdt(ip)=(press0(ip)-press1(ip))/stepseconds
          do iz=0,npro
            dthdt(iz,ip)=(thpro0(iz,ip)-thpro1(iz,ip))/mrefst
            dthdzdt(iz,ip)=(dthdz(iz,ip)-dthdz1(iz,ip))/mrefst
            dusdt(iz,ip)=(uspro0(iz,ip)-uspro1(iz,ip))/mrefst
            dusdzdt(iz,ip)=(dusdz(iz,ip)-dusdz1(iz,ip))/mrefst
            dvsdt(iz,ip)=(vspro0(iz,ip)-vspro1(iz,ip))/mrefst
            dvsdzdt(iz,ip)=(dvsdz(iz,ip)-dvsdz1(iz,ip))/mrefst
            dqvsdt(iz,ip)=(qvspro0(iz,ip)-qvspro1(iz,ip))/mrefst
            dqvsdzdt(iz,ip)=(dqvsdz(iz,ip)-dqvsdz1(iz,ip))/mrefst
          enddo
        enddo

! verify mass balance in the reference state profiles
! computations at time zero (of the update cycle)

        if(divrefstate.and. nprof.gt.1) then

          dppref=0.
          dpmeddt=0.
          do ip=1,nprof
            do ix=0,nx1
            do iy=0,ny1
              dppref(ix,iy)=dppref(ix,iy)+pweight(ix,iy,ip)*dpressdt(ip)
            enddo
            enddo
            do ix=2,nx
            do iy=2,ny
              dpmeddt=dpmeddt+dppref(ix,iy)
            enddo
            enddo
          enddo

          if(.not.allocated(ubxpro1)) allocate(ubxpro1(2,2:ny,0:npro))
          if(.not.allocated(ubxpro2)) allocate(ubxpro2(2,2:ny,0:npro))
          if(.not.allocated(dubxprodt))
     :      allocate(dubxprodt(2,2:ny,0:npro))
          if(.not.allocated(dubxprodz1))
     :      allocate(dubxprodz1(2,2:ny,0:npro))
          if(.not.allocated(dubxprodz2))
     :      allocate(dubxprodz2(2,2:ny,0:npro))
          if(.not.allocated(dubxprodzdt))
     :      allocate(dubxprodzdt(2,2:ny,0:npro))
          if(.not.allocated(rhox)) allocate(rhox(2,2:ny,0:npro))
          if(.not.allocated(xflop)) allocate(xflop(2,2:ny,0:npro))

          if(.not.allocated(vbypro1)) allocate(vbypro1(2:nx,2,0:npro))
          if(.not.allocated(vbypro2)) allocate(vbypro2(2:nx,2,0:npro))
          if(.not.allocated(dvbyprodt))
     :      allocate(dvbyprodt(2:nx,2,0:npro))
          if(.not.allocated(dvbyprodz1))
     :      allocate(dvbyprodz1(2:nx,2,0:npro))
          if(.not.allocated(dvbyprodz2))
     :      allocate(dvbyprodz2(2:nx,2,0:npro))
          if(.not.allocated(dvbyprodzdt))
     :      allocate(dvbyprodzdt(2:nx,2,0:npro))
          if(.not.allocated(rhoy)) allocate(rhoy(2:nx,2,0:npro))
          if(.not.allocated(yflop)) allocate(yflop(2:nx,2,0:npro))


          ubxpro1=0.
          rhox=0.
          vbypro1=0.
          rhoy=0.

          do ip=1,nprof
            do iz=0,npro
            do iy=2,ny
              ubxpro1(1,iy,iz)=ubxpro1(1,iy,iz)
     :          +pweight(1,iy,ip)*(uspro1(iz,ip)+dusdz1(iz,ip)*zpro(iz))
              ubxpro1(2,iy,iz)=ubxpro1(2,iy,iz)
     :          +pweight(nx,iy,ip)*(uspro1(iz,ip)
     :          +dusdz1(iz,ip)*zpro(iz))

              rhox(1,iy,iz)=rhox(1,iy,iz)
     :          +pweight(1,iy,ip)*rhopro1(iz,ip)
              rhox(2,iy,iz)=rhox(2,iy,iz)
     :          +pweight(nx,iy,ip)*rhopro1(iz,ip)
            enddo
            enddo

            do iz=0,npro
            do ix=2,nx
              vbypro1(ix,1,iz)=vbypro1(ix,1,iz)
     :          +pweight(ix,1,ip)*(vspro1(iz,ip)+dvsdz1(iz,ip)*zpro(iz))
              vbypro1(ix,2,iz)=vbypro1(ix,2,iz)
     :          +pweight(ix,ny,ip)*(vspro1(iz,ip)
     :          +dvsdz1(iz,ip)*zpro(iz))

              rhoy(ix,1,iz)=rhoy(ix,1,iz)
     :          +pweight(ix,1,ip)*rhopro1(iz,ip)
              rhoy(ix,2,iz)=rhoy(ix,2,iz)
     :          +pweight(ix,ny,ip)*rhopro1(iz,ip)
            enddo
            enddo
          enddo

          flowinx=0.
          flowoutx=0.
          flowiny=0.
          flowouty=0.

          do iy=2,ny
          do iz=0,npro
            xflop(1,iy,iz)=ubxpro1(1,iy,iz)*rhox(1,iy,iz)
            xflop(2,iy,iz)=ubxpro1(2,iy,iz)*rhox(2,iy,iz)
          enddo
          enddo

          do iy=2,ny
          do iz=0,npro
            flowinx=flowinx+max(0.d0,xflop(1,iy,iz))
     :         -min(0.d0,xflop(2,iy,iz))
            flowoutx=flowoutx-min(0.d0,xflop(1,iy,iz))
     :         +max(0.d0,xflop(2,iy,iz))
          enddo
          enddo

          do iz=0,npro
          do ix=2,nx
            yflop(ix,1,iz)=vbypro1(ix,1,iz)*rhoy(ix,1,iz)
            yflop(ix,2,iz)=vbypro1(ix,2,iz)*rhoy(ix,2,iz)
          enddo
          enddo

          do iz=0,npro
          do ix=2,nx
            flowiny=flowiny+max(0.d0,yflop(ix,1,iz))
     :        -min(0.d0,yflop(ix,2,iz))
            flowouty=flowouty-min(0.d0,yflop(ix,1,iz))
     :        +max(0.d0,yflop(ix,2,iz))
          enddo
          enddo

          flowin=flowinx+flowiny
          flowout=flowoutx+flowouty
          area=(nx-1)*dx*(ny-1)*dy

          residual=flowin-flowout-dpmeddt*area
          totalcf=flowin+flowout

          if(residual.gt.1.e-6) then
            do iy=2,ny
            do iz=0,npro
              ubxpro1(1,iy,iz)=ubxpro1(1,iy,iz)
     :          -residual*abs(xflop(1,iy,iz))/(totalcf*rhox(1,iy,iz))
              ubxpro1(2,iy,iz)=ubxpro1(2,iy,iz)
     :          -residual*abs(xflop(2,iy,iz))/(totalcf*rhox(2,iy,iz))
            enddo
            enddo

            do iz=0,npro
            do ix=2,nx
              vbypro1(ix,1,iz)=vbypro1(ix,1,iz)
     :          -residual*abs(yflop(ix,1,iz))/(totalcf*rhoy(ix,1,iz))
              vbypro1(ix,2,iz)=vbypro1(ix,2,iz)
     :          -residual*abs(yflop(ix,2,iz))/(totalcf*rhoy(ix,2,iz))
            enddo
            enddo
          endif


! veri  fy mass balance in the reference state profiles
! comp  utations at time max (of the update cycle)

          ubxpro2=0.
          vbypro2=0.
          rhox=0.
          rhoy=0.

          do ip=1,nprof

            do iz=0,npro
            do iy=2,ny
              ubxpro2(1,iy,iz)=ubxpro2(1,iy,iz)
     :          +pweight(1,iy,ip)*(uspro0(iz,ip)+dusdz(iz,ip)*zpro(iz))
              ubxpro2(2,iy,iz)=ubxpro2(2,iy,iz)
     :          +pweight(nx,iy,ip)*(uspro0(iz,ip)+dusdz(iz,ip)*zpro(iz))

              rhox(1,iy,iz)=rhox(1,iy,iz)
     :          +pweight(1,iy,ip)*rhopro(iz,ip)
              rhox(2,iy,iz)=rhox(2,iy,iz)
     :          +pweight(nx,iy,ip)*rhopro(iz,ip)
            enddo
            enddo

            do iz=0,npro
            do ix=2,nx
              vbypro2(ix,1,iz)=vbypro2(ix,1,iz)
     :          +pweight(ix,1,ip)*(vspro0(iz,ip)+dvsdz(iz,ip)*zpro(iz))
              vbypro2(ix,2,iz)=vbypro2(ix,2,iz)
     :          +pweight(ix,ny,ip)*(vspro0(iz,ip)+dvsdz(iz,ip)*zpro(iz))

              rhoy(ix,1,iz)=rhoy(ix,1,iz)
     :          +pweight(ix,1,ip)*rhopro(iz,ip)
              rhoy(ix,2,iz)=rhoy(ix,2,iz)
     :          +pweight(ix,ny,ip)*rhopro(iz,ip)
            enddo
            enddo
          enddo

          flowinx=0.
          flowoutx=0.
          flowiny=0.
          flowouty=0.

          do iy=2,ny
          do iz=0,npro
            xflop(1,iy,iz)=ubxpro2(1,iy,iz)*rhox(1,iy,iz)
            xflop(2,iy,iz)=ubxpro2(2,iy,iz)*rhox(2,iy,iz)
          enddo
          enddo

          do iy=2,ny
          do iz=0,npro
            flowinx=flowinx+max(0.,xflop(1,iy,iz))
     :         -min(0.,xflop(2,iy,iz))
            flowoutx=flowoutx-min(0.,xflop(1,iy,iz))
     :         +max(0.,xflop(2,iy,iz))
          enddo
          enddo

          do iz=0,npro
          do ix=2,nx
            yflop(ix,1,iz)=vbypro2(ix,1,iz)*rhoy(ix,1,iz)
            yflop(ix,2,iz)=vbypro2(ix,2,iz)*rhoy(ix,2,iz)
          enddo
          enddo

          do iz=0,npro
          do ix=2,nx
            flowiny=flowiny+max(0.,yflop(ix,1,iz))
     :        -min(0.,yflop(ix,2,iz))
            flowouty=flowouty-min(0.,yflop(ix,1,iz))
     :        +max(0.,yflop(ix,2,iz))
          enddo
          enddo

          flowin=flowinx+flowiny
          flowout=flowoutx+flowouty
          area=(nx-1)*dx*(ny-1)*dy

          dppref=0.
          psufref=0.
          dpmeddt=0.

c note  : 1/g cancels

          residual=flowin-flowout-dpmeddt*area
          totalcf=flowin+flowout

c         write(*,*) 'residual:',residual,flowin,flowout,dpmeddt*area

          if(residual.gt.1.e-6) then
            do iy=2,ny
            do iz=0,npro
              ubxpro2(1,iy,iz)=ubxpro2(1,iy,iz)
     :          -residual*abs(xflop(1,iy,iz))/(totalcf*rhox(1,iy,iz))
              ubxpro2(2,iy,iz)=ubxpro2(2,iy,iz)
     :          -residual*abs(xflop(2,iy,iz))/(totalcf*rhox(2,iy,iz))
            enddo
            enddo

            do ix=2,nx
            do iz=0,npro
              vbypro2(ix,1,iz)=vbypro2(ix,1,iz)
     :          -residual*abs(yflop(ix,1,iz))/(totalcf*rhoy(ix,1,iz))
              vbypro2(ix,2,iz)=vbypro2(ix,2,iz)
     :          -residual*abs(yflop(ix,2,iz))/(totalcf*rhoy(ix,2,iz))
            enddo
            enddo
          endif

          do ii=1,2
          do iy=2,ny
          do iz=0,npro
            dubxprodt(1,iy,iz)=(ubxpro2(1,iy,iz)-ubxpro1(1,iy,iz))
     :        /mrefst
          enddo
          enddo
          enddo

          do ii=1,2
          do iy=2,ny
          do iz=0,npro-1
            dubxprodz1(ii,iy,iz)=(ubxpro1(ii,iy,iz+1)-ubxpro1(ii,iy,iz))
     :        /dzpro
            dubxprodz2(ii,iy,iz)=(ubxpro2(ii,iy,iz+1)-ubxpro2(ii,iy,iz))
     :        /dzpro
          enddo
          enddo
          enddo

          do ii=1,2
          do iy=2,ny
          do iz=0,npro-1
            dubxprodzdt(ii,iy,iz)=(dubxprodz2(ii,iy,iz)
     :        -dubxprodz1(ii,iy,iz))/mrefst
          enddo
          enddo
          enddo

          do ii=1,2
          do ix=2,nx
          do iz=0,npro
            dvbyprodt(ix,ii,iz)=(vbypro2(ix,ii,iz)-vbypro1(ix,ii,iz))
     :        /mrefst
          enddo
          enddo
          enddo

          do ii=1,2
          do ix=2,nx
          do iz=0,npro-1
            dvbyprodz1(ix,ii,iz)=(vbypro1(ix,ii,iz+1)
     :        -vbypro1(ix,ii,iz))/mrefst
            dvbyprodz2(ix,ii,iz)=(vbypro2(ix,ii,iz+1)
     :        -vbypro2(ix,ii,iz))/mrefst
          enddo
          enddo
          enddo

          do ii=1,2
          do ix=2,nx
          do iz=0,npro-1
            dvbyprodzdt(ix,ii,iz)=(dvbyprodz2(ix,ii,iz)
     :        -dvbyprodz1(ix,ii,iz))/dzpro
          enddo
          enddo
          enddo

          if(.not.allocated(dweightx)) then
            allocate(dweightx(0:nx+1,0:ny+1,2))
            allocate(dweighty(0:nx+1,0:ny+1,2))
            do ix=0,nx+1
              distx2=(ix-nx)*dx
              distx1=(ix-2)*dx
              do iy=0,ny+1
                disty2=(iy-ny)*dy
                disty1=(iy-2)*dy
                distinv11=1./sqrt(distx1**2+disty1**2)
                distinv12=1./sqrt(distx1**2+disty2**2)
                distinv21=1./sqrt(distx2**2+disty1**2)
                distinv22=1./sqrt(distx2**2+disty2**2)
                sumdisti=distinv11+distinv12+distinv21+distinv22
                dweightx(ix,iy,1)=distinv11/sumdisti
                dweightx(ix,iy,2)=distinv21/sumdisti
                dweighty(ix,iy,1)=distinv12/sumdisti
                dweighty(ix,iy,2)=distinv22/sumdisti
              enddo
            enddo
          endif
        endif

! end_divrefstate

        if(verbose.ge.2) then
          do ip=1,nprof
            write(*,*) 'New profile',ip,dpressdt(ip)
            do iz=1,npro
            write(*,'(i5,8e15.7)')
     :        iz,thpro0(iz,ip)
     :        ,uspro0(iz,ip),vspro0(iz,ip)
     :        ,qvspro0(iz,ip)
            enddo
          enddo
        endif

c compute tendencies

        do ip=1,nprof
          do iz=0,npro
            dusdt(iz,ip)=(uspro0(iz,ip)-uspro1(iz,ip))/mrefst
            dusdzdt(iz,ip)=(dusdz(iz,ip)-dusdz1(iz,ip))/mrefst
            dvsdt(iz,ip)=(vspro0(iz,ip)-vspro1(iz,ip))/mrefst
            dvsdzdt(iz,ip)=(dvsdz(iz,ip)-dvsdz1(iz,ip))/mrefst
          enddo
        enddo

c recover saved profiles

        thpro0=thpro1
        uspro0=uspro1
        vspro0=vspro1
        qvspro0=qvspro1
        press0=press1

        dthdz=dthdz1
        dusdz=dusdz1
        dvsdz=dvsdz1
        dqvsdz=dqvsdz1

        if(verbose.ge.2) then
          do ip=1,nprof
            write(*,'(a30,i5,2e15.7)') 'Uprefst Next current profile'
     :        ,ip,press0(ip),dpressdt(ip)
            do iz=1,npro
            write(*,'(i5,8e15.7)')
     :        iz,thpro0(iz,ip),dthdt(iz,ip)
     :        ,uspro0(iz,ip),dusdt(iz,ip),vspro0(iz,ip),dvsdt(iz,ip)
     :        ,qvspro0(iz,ip),dqvsdt(iz,ip)
            enddo
          enddo
        endif

      endif


      do ip=1,nprof
        press0(ip)=press0(ip)+dpressdt(ip)*dt
        do iz=0,npro
          thpro0(iz,ip)=thpro0(iz,ip)+dthdt(iz,ip)
          dthdz(iz,ip)=dthdz(iz,ip)+dthdzdt(iz,ip)
          uspro0(iz,ip)=uspro0(iz,ip)+dusdt(iz,ip)
          dusdz(iz,ip)=dusdz(iz,ip)+dusdzdt(iz,ip)
          vspro0(iz,ip)=vspro0(iz,ip)+dvsdt(iz,ip)
          dvsdz(iz,ip)=dvsdz(iz,ip)+dvsdzdt(iz,ip)
          qvspro0(iz,ip)=qvspro0(iz,ip)+dqvsdt(iz,ip)
          dqvsdz(iz,ip)=dqvsdz(iz,ip)+dqvsdzdt(iz,ip)
        enddo
      enddo

      if(verbose.ge.3) then
        do ip=1,nprof
          write(*,'(a20,i6,e15.7)') 'New profile',ip,press0(ip)
          do iz=1,npro
          write(*,'(i5,4f10.3)')
     :      iz,thpro0(iz,ip),uspro0(iz,ip),vspro0(iz,ip),qvspro0(iz,ip)
          enddo
        enddo
      endif

      return
8888  continue
      write(0,*) 'nh3dFatalError=ERROR in uprefst (end of file)'
      stop
      end
